function [J0,g_norm,sol,delta_sol,saddle,u_max]=minimax(n,u,delta_u,solution) 
global c a_t b p e t n_t n_p x_t y_t ar nth0 pwr;

epsilon=input('Starting Minimax Method. Input its accuracy = ');
%      s_max=input('Input maximum stepsize s_max=');
       s_max=0.01;
       it_max=input('Input Maximum Iteration Number it_max=');

if n==1    %MI=1
hom=input('Is the nonlinear term homogeneous? 1-yes,0-no=');
end

disp('     It_m     g_norm    err_max   ener gy    u_max at       (x,y)           s');

r=zeros(n,1);

if (n==1 & hom==1)

top=ar*((-1/c*delta_u(:,n)+a_t.*u(:,n)).*u(:,n));
bot=ar*u(:,n).^(pwr+1);
r(n)=(top/bot)^(1/(pwr-1));
r0=r;
J0=(0.5*r(n)^2*top-r(n)^(pwr+1)/(pwr+1)*bot);

else

r(n)=1.0;
r(n)=fminsearch('energy',r(n),[],u(:,n),delta_u(:,n),1);
%*************************************************************
options=optimset('GradObj','on','Display','off','TolX',1e-6);
[r0,J0]=fminunc('energy',r,options,u,delta_u,n);
%**************************************************************
%[r0,J0]=fminsearch('energy',r,[],u,delta_u,n);
J0=-J0;

end

[g_t,delta_g_t,g_p,g_norm,res,pu,delta_pu]=gradient(r0,u,delta_u,n);
%pu=p(u) peak slection, res=residual of PDE, g_t,g_p=negative gradient J,

J=J0;
it_n=1;
w=u;
delta_w=delta_u;

while g_norm>epsilon
      s=s_max;
%      while J >= J0
       while J-J0>= -0.5*r(n)*s*g_norm^2     %Check stepsize rule
    	w(:,n)=u(:,n)+s*g_t';
	    delta_w(:,n)=delta_u(:,n)+s*delta_g_t';

if (n==1 & hom==1)

top=ar*((-1/c*delta_w(:,n)+a_t.*w(:,n)).*w(:,n));
bot=ar*w(:,n).^(pwr+1);
r(n)=(top/bot)^(1/(pwr-1));

J=(0.5*r(n)^2*top-r(n)^(pwr+1)/(pwr+1)*bot);

else
%*************************************************************
options=optimset('GradObj','on','Display','off','TolX',1e-6);
[r,J]=fminunc('energy',r0,options,w,delta_w,n);
%**************************************************************
%        [r,J]=fminsearch('energy',r0,[],w,delta_w,n);
        J=-J;
end

        s=s/2;
        if (s <= 0.00001 & J < J0)           %relax stepsize rule
           break;
        end
      end
      u(:,n)=w(:,n);
      delta_u(:,n)=delta_w(:,n);
      solution(:,n)=solution(:,n)+2*s*g_p;
      J0=J;
      r0=r;
      [g_t,delta_g_t,g_p,g_norm,res,pu,delta_pu]=gradient(r0,u,delta_u,n);
      err_max=max(abs(res));
      [u_max, II]=max(abs(pu));
      disp([it_n g_norm err_max J0 u_max x_t(II) y_t(II) s]);

%      pu_p=zeros(n_p,1);
%      for i=1:n
%       pu_p=pu_p+r0(i)*solution(:,i);
%      end
      pu_p = solution*r0;

      u_max=-min(pu_p);
      u_min=max(pu_p);
      set(gca,'ZLim',[u_min u_max]);
      %pdesurf(p,t,pu_p); drawnow;
      pdeplot(p,[],t,'XYData',pu_p,'ZData',-pu_p,'ColorBar','off'); drawnow;
      it_n=it_n+1;
      if it_n-1 > it_max
      ans=input('It_n > It_max, want more Minimax iteration? 1-yes, 0-no=');
       if ans==1
          it_max=it_max+input('Input the number of MORE Minimax iterations=');
       else
          break;
       end
      end
      s_max=3*s;
  end
%sol=zeros(n_t,1);
%delta_sol=zeros(n_t,1);
%saddle=zeros(n_p,1);
%for i=1:n
%  sol=sol+r0(i)*u(:,i);
%  delta_sol=delta_sol+r0(i)*delta_u(:,i);
%  saddle=saddle+r0(i)*solution(:,i);
%end
sol = u*r0;
delta_sol = delta_u*r0;
saddle = solution*r0;

